
p <- matrix(c(4, 3, 2, 1, 3, 4, 3, 2, 2, 3, 4, 3, 1, 2, 3, 4) / 40, 4, 4, byrow =
              TRUE)


print(p)

print(dim(p))

print(sum(p))

dimnames(p)[[1]] <- 1:4
dimnames(p)[[2]] <- 1:4

p[, 1]
p[2, ]

niter <- 1000

gibbs_discrete <- function(p, i = 1, iter = niter) {
  x <- matrix(0, iter, 2)
  nX <- dim(p)[1]
  nY <- dim(p)[2]
  
  for(k in 1:iter){
    j <- sample(1:nY, 1, prob=p[i,])
    i <- sample(1:nX, 1, prob=p[,j])
    
    x[k,] <- c(i,j)
  }
  
  x
}

sp <- data.frame(gibbs_discrete(p))

names(sp) <- c('X','Y')
print(p)

table(sp)/niter

# prop.table(sp)



dbinom(10,20, 0.5)

pbinom(15, 20, 0.5)


rbinom(1,size=20, prob=0.5)

